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Image magnification by least squares 
surfaces 


A.M. Esmaili Zaini, G.B. Loghmani*, and A.M. Latif 


Abstract 


Image magnification is one of the current issues of image processing in 
which keeping the quality and structure of images is the main concern. In im- 
age magnification, it is necessary to insert information in extra pixels. Adding 
information to an image should be compatible with the image structure with- 
out making artificial blocks. In this research, extra pixels are estimated using 
the surface of least squares, and all the pixels are reviewed according to the 
suggested edge-improving algorithm. The suggested method keeps the edges 
and minimizes the magnified image opacity and the artificial blocks. Numer- 
ical results are presented by using PSNR and SSIM fidelity measures and 
compared to some other methods. The average PSNR of the original image 
and image zooming is 32.79 which it shows that image zooming is very similar 
to the original image. Experimental results show that the proposed method 
has a better performance than others and provides good image quality. 
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1 Introduction 


Image magnification is, indeed, image resolution increasing to achieve a 
higher- quality image. It plays an important role in image processing and 
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machine vision. Image magnification has various applications in electronic 
publishing industries, digital cameras, medical imaging, images on the web, 
license plate recognition, and face recognition systems in law functions. 

Recent studies indicate that the main emphasis is on the visual quality 
of images in many applications. Resolution edges and lack of blur and addi- 
tional artifacts are two important factors in image quality. Many enlarging 
image algorithms use interpolation methods. Interpolation means to find a 
set of unknown pixel values in a set of known pixel values in the picture. 
In magnification, the following few basic parameters affect the image qual- 
ity [11]: 

1- An efficient magnification method should preserve the edges and bound- 
aries. 

2- The method should not produce undesirable constant piecewise or blocks 
of other regions. 

3- The magnification method should be computationally efficient but not too 
dependent on the internal parameters of the image. 

Traditional technologies in image magnification use linear interpola- 
tion for high-resolution samples. Pixel replication, bilinear interpolation, 
quadratic interpolation, bi-cubic interpolation, and spline interpolation are 
some of these methods [6, 8,17]. These methods tend to smooth edges or 
produce blocky interpolated images with staircase edges. Therefore, the out- 
put of these methods produces blurred images. This indicates the inability 
of linear technologies to transfer new information to an image. 

One of the features of the two linear interpolation methods is that, in 
the magnification ratio, artificial blocks and visual effects are undesirable, 
but edges are preserved at an acceptable level. In bi-cubic interpolation 
method with a high zoom ratio, artificial blocks and undesirable visual effects 
are lower, and the edges are preserved. Although determining the image 
quality is not easy in this method, for image magnification, the quadratic 
interpolation is better than bilinear interpolation, bi-cubic interpolation is 
better than quadratic interpolation, and spline interpolation is better than 
bi-cubic interpolation [16]. 

In recent years, nonlinear interpolation methods have been used to reform 
linear methods for improving image quality and solving the blur problem. 
The change of non-linear methods depends on the interpolation method. This 
means that the performance of these methods with sharp edges is different 
from their performance on soft tissues, while linear methods with all the 
pixels act in the same way [1,3,7]. In order to maintain the image quality 
and edges in non-linear methods, estimation methods of a subset of edge 
pixels [1], resampling and optimization parameters, implicit interpolation, 
and straight edges have been used [4, 14, 22]. 

In [22], a statistical method that tunes interpolation coefficients according 
to local edge structures is proposed. The technique of [4] uses a modified 
bilinear method, where the interpolation error theorem in an edge-adaptive 
fashion. 
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In [21], an edge-directed nonlinear interpolation technique is presented 
through directional filtering and data fusion. This algorithm interpolates a 
missing pixel in multiple directions, and then fuses the directional interpola- 
tion results by linear minimum mean square-error estimation. 


In [2], the local adaptive magnification (LAZ) method uses the discon- 
tinuities information or sharp brightness changes, finds edges in different 
directions by taking two threshold values, and estimates the unknown pixels 
with respect to the edge pixels. 


In [20], the artificial neural networks methods uses for doing zoom oper- 
ations on digital images. 


In [10], presents a nonlinear image interpolation algorithm that is based on 
the moving least squares technique. This methods ability to image zooming 
and preserving edge features. 


In this paper, another least squares surface method is suggested, in which 
necessary pixels are obtained by calculating their coefficients using the least 
squares theory and edge-directed algorithms. 


This paper continues as follows. In the second part, quadratic surfaces 
and the theory of least squares will be discussed. In the third part, the 
least square planes, suggested algorithms, and evaluation parameters will be 
proposed. The fourth section compares the results of implementation by 
other methods. In the last section, conclusions and recommendations will be 
presented. 


2 Quadratic surfaces and least squares theory 


The purpose of a quadratic surface in R® space composed of x, y and z is 
formed as equation (1), in which none of the coefficients A, B, C, D, E and 
F is zero. This equation called quadratic surfaces is an extension of a cone 
in R? space. 


Az? + By? + Cz? 4+ Day4+ Exz4+ Fyz+Gr+Hy+Iz+J=0. (1) 


Without losing any generality, the coefficients of product terms can be 
zero by a three-dimensional rotation, and convert equation (1) into a conven- 
tional or canonical form by a transfer, whose diagram is recognizable. The 
canonical equation for a quadratic surface is as in (2), the most famous sec- 
tions of which are elliptical, hyperbolic, parabolic, and conic sections. 

For example, the overall shape of an elliptical equation and parabolic surface 
are introduced according to equations (3) and (4), respectively: 
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Az? + By? +Cz?+Gr+ Hy+Iz+J=0 (2) 
Az? + By? +C2?4+Gze+Hyt+Iz=1 (3) 
Az? + By? ++G2 + Hy = Iz (4) 


2.1 Least Squares Theory 


Suppose f(x) is a continuous function on the [a, 6] interval and p(x) is a 
polynomial of the maximum degree n that is defined according to equation 
(5), in which c; constants for i = 0,1,2,...,n are real numbers. 


p(a) = os Cu". (5) 


The problem of approximating function f(x) by polynomial p(x) can be 
considered by finding c; constants for i = 0,1,2,...,n to minimize || f(a) — 
p(z)||. The purpose of |].|| is Euclidean norms. To calculate the coefficients 
of c;, Theorem 1 is used. 


Theorem 1. Let f (21, 2%2,...,%n) be a function of n variables which has 
a local extremum at the point (a,,a2,...,@,). Then either f is non- dif- 
ferentiable at (a1, 42,...,@n) , or it is differentiable at (a1, a2,...,Qn) and 


Vf (ai, @2,---,@n) =0 or Fay at 2s Gn) = 0, CHORD A 


If a set of data are in the form of {(a;, y:) |i = 1,2,...,.N} and the suitable 
mathematical model for data relationship is a line, then the general problem is 
finding the best least squares line. A mathematical model for the relationship 
of the data is the polynomial of degree n in general for which n < N., Then 
the problem is to find the least squares polynomial. Therefore, if the discrete 
data based on Euclidean norms are considered as in expression (6) according 
to Theorem 1, equation (7) should be established. 


N 
E= (ys — (co teyay +... + en24"))? (6) 
i=l 
OE 
— —0 7=0,1,2,...,7. 7 
Oc; F J 2 9 F 3 Pa (7) 


After calculating the partial derivatives and simplifying them, the linear sys- 
tem of (n+ 1) equations called normal form is obtained as in (8) for distinct 
x;, fori =1,2,...,N. This linear system of equations has a unique solution. 
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N N N N 
cod, ee ted att... tend) 2? = >> yx? 
i=l i=l j= i= 


N N ‘y , 
cod. tt tad att... tend ct = ¥Y yal 
i=l i=l i=l i=1 (8) 
ae oer "Reva ch 
com ep tery ap +...+en >) 27" = do yah 
i=1 i=1 i=1 i=1 


If f(z, y,z) = 0 is a surface equation in the R® space, z depends on x and 
y. In addition, {(ai, yi, zi) |i = 1,2,...,N} is a set of discrete data obtained 
from the process of data collection, non of which lie on a known surface, 
either elliptical, or parabolic, or plane. Then according to 1 and the sum 
squares error, appropriate coefficients can be calculated. 


Suppose z = Ax + By+C is a surface equation in the R® space for 
{(xi, Yi, 2%) |i =1,2,...,N} points. To obtain the surface of least squares 
Theorem 1 and error (9), the coefficients of A, B,C are calculated after solv- 
ing normal equations (10). 


E= 


4 


(Ac, + By, + C— %) ? (9) 


N 
=1 


N N N 
Ad a7 + BY ey + CyS ap = >) aa; 
i=1 i= 1 =I 


N Me N v 
Ad tit Bw + Ch w= Li wie (10) 


N N N N 
Adi «t,+ BO y+CO1= Su. 
i=l i=l i=l i=l 


If Ax? + By? + Cx + Dy = 1-— 2? is an elliptical surface in the space R* and 
{(xi, yi, 21) |i =1,2,...,N} points are provided, calculation of the elliptical 
least-squares surface can be done according to Theorem 1 and deviation (11), 
with coefficients A, B,C, D achieved after solving linear device (12). 


E= 


a 


(Ax? + By? + Ca;+ Dy; + a - 1) 2 (11) 


N 
= 
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N N N N N 
A>) a$ + BY 227y? +CY) 23 + DDS yx? = YS 2? (1 — 2?) 
= i=l 


w=1 w=1 i=1 


2 


. 29,2. a 4 of 2 = af 2 2 
eae = i= 
N 


i=l i=l 


N N N N 
Ay a3 + BY eyZ+Cd a2 + DY wa, = ¥ x; (1-2?) 
i i=l i=l i=l 


i= i=1 


Rs N N Ms N : 
Ay wut BY y+ Cd sy, + Dd yy = 2 Yi (127): 


i= 


(12) 


° 


Similarly, if Ax? + By? + Ca + Dy = z is parabolic surface in the space 
R® and {(2j, yi, zi) |i = 1,2,...,N} points are provided, calculation of the 
parabolic least-squares surface can be done according to Theorem 1 and de- 
viation (13), with coefficients A, B,C, D achieved after solving linear device 
(14). 


E= 


a 


(Ax? + By? + Ca; + Dy, — wy (13) 


N 
Al 


N N N N N 
A>) ap + BY) aly} + Cyl a3 + DY ya? = Do 272, 
i=l i=l i=l i=l io 


N N N N 
AX 22y?+ BY yf +CD ay? + DY B= YT yee; 
i=l i=l i=l i=l i=l 
N N N N N 
AY, BY tye PCy a} PD a= a 
i=l i=l i=l i=l i=l 


W Me N N. ON 
AD) tiyt BY yg t+ Cyl ay, + DY yf = Ye 
i=l =1 i=1 ral i=l 


(14) 


uv 


Here, the linear equations system can be obtained for the quadratic surface. 


3 Least Squares Surfaces and Image Magnification 


In the image magnification, some new pixels are placed in the original pixels 
of the image. The purpose of magnification is to determine new pixels, which 
are determined based on their neighbouring pixels. There are two types of 
neighbourhood in two-dimensional images, which are known as the 4 and 8 
cell. Neighbourhoods are shown in Figure 1. 


In most of proposed methods, magnification rate is tried to be squared, 
or to a power of two, but it should be noted that the algorithm presented 
here is used for every magnification rate. 

The first stage is the simplest one and requires expanding the source n x n 
pixels image onto a regular grid of size (2n — 1) x (2n —1). More precisely, if 
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1 2 3 


4" neighborhoods 8" neighborhoods 


Figure 1: Types of neighbourhood. 


S(i,7) denotes the pixels in the i” row and j‘” column of the source image 
and Z(1,k) denotes the magnified image pixel in the I'” row and k*” column 
in the zoomed picture, the f function [5] puts the original image pixel in the 
interlaced places of the new image: 


f:S7Z 
f (S(é,7)) = Z (2 -— 1,27 - 1) sy 9 Sy Qy ccd pit 


The result is shown in Figure 2. The original image pixel is indicated by 
symbol e, and the other pixels which must be estimated have been identified 
by 11,10 and 01. To estimate the pixel values, the multi-stage least squares 
method is used along with a review of the surfaces. 

In this algorithm, a new pixel value is estimated with its four neigh- 
bourhoods in the corner (pixel 11) by the proposed method. Since, all the 
neighbourhoods on the central pixels are equal in term of distance from the 
central pixel, all of them are attributed the same weight. The results are 
shown in Figure 3, and M is replaced by the estimated pixels. In the next 
stage, the pixels indicated by 1 are estimated by their four neighbourhoods. 
These pixels have two main pixels in the left and the right and two esti- 
mated pixels up and down the neighbourhoods. In estimation, pixels that 
are indicated with number 01 have the weight of 1, and the estimated pixels 
assigned with the weight of 0.5. Pixels with number 10 are estimated in the 
same manner. The results are shown in Figure (4), in which the symbols 
and are replaced by the estimated pixels. 


3.1 Least squares surface algorithm 


In this algorithm, a point (e) refers to the pixel coordinates of the image 
whose brightness and number of rows and columns are considered as a point 
in a three-dimensional space. 
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Figure 4: Estimation of neighbourhood pixels. 


To estimate the required pixels in image magnification based on four 
neighbourhoods, initially by selection the four points and solving the linear 
least-squares surface equation, the desired pixel is estimated. It should be 
noted that, two adjacent pixels are used in the rows and columns of pixels to 
approximate the first and last of the average. 

The proposed algorithm steps are as follows: 

Step 1. Establish a linear system according to the chosen least-square sur- 
face. 

Step 2. Solve linear equations and calculating surface coefficients. 
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Step 3. Estimate the desired pixel by replacing the number of rows and 
columns on the selective surface. 

Step 4. Repeat steps 1 to 3 in three stages according to the zooming algo- 
rithm and estimating the required pixel 

Step 5. Review all the estimated pixels based on the conducted sub-algorithm 
of the edge (3-2). 


3.2 The modified sub-algorithm of the image edges 


After estimating all the required image pixels, all estimated pixels are re- 
viewed to improve the edges (Figure.5) based on the following conditions.A, B, 
C, and d are the main pixels, and X,Y2,Y1, 71, and Z2 are the estimated 
pixels. The advantage of the sub-algorithm is the lack of any provinces for 
edges. The proposed algorithm steps are as follows: 

Step 1. If |A—D| > |B—C|, then X «+ 2£° | In fact, the edge is in the 
northeast-southwest direction. 

Step 2. If |A—D| <|B-—C|, then X + 44". In fact, the edge is in the 
northwest-southeast direction. 

Step 3. If (A— D)(B-—C) > 0, then Y1 + 444 and Y2+ SS | In fact, 
the edge is in the north-south direction. 

Step 4. If (A— D)(B-—C) <0, then Z1 + “$£and Z2 + 24 | In fact, 
the edge is in the west-east direction. 


a. | ee 
ZL |*% |22 
Ce |*2 | 


Figure 5: Design for modifying the image edge. 


3.3 Assessment criteria 


To compare two images in terms of their similarity, PSNR and SSIM criteria 
are used. The PSNR value is bigger and the SSIM is closer to 1, the images 
have more conformity [18,19]. The PSNR criterion is calculated based on 
equation (15) and (16): 
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MAX? 
PSNR = 10 Logio (a) (dB) (15) 
1 m—-l1n-1 
MSE=— 0 DU GA) -—K GA (16) 
i=0 j=0 


where J (i, 7) and K (%,7) are the main image pixels and the estimated image, 
respectively, and MAX, is the maximum image pixles. 

SSIM criterion, which includes the image structural elements and measures 
the structural quality of the image as well as the similarity between the two 
images, has a value between 0 and 1. This means that 1 is the highest and 0 
is the lowest similarity, which is calculated by equation (17). 


(pte by + C1) (27 2y + €2) 


IM =F 
SS (x,y) (u2 am pe + c1) (0? + of + c2) 


(17) 


where yp , 07 and oxy are the average, variance, and covariance of the pixels 
in the image, respectively, and cj and cp are two unknown constants that 
prevent the fraction denominator from being zero. 


4 Simulation results 


To evaluate the proposed method, first, a digital image is considered as an 
original image, and then its size is reduced by removing alternative rows of 
columns by half, then, it is tried to make image size double using the pro- 
posed method and the other methods. As expected, the closer the degree 
of similarity is to the original image size, the algorithm has a better perfor- 
mance. 


Implementation of the results of the suggested algorithm done in MAT- 
LAB. The results are presented in Tables 1 and 2. In Table 1, the PSNR 
criterion is used for zooming, and ten different and standard images [15] are 
applied according to Figure 6. Their dimensions all are 512 x 512. The results 
of the suggested least square plane (LSP) and least square ellipsoid (LSE) 
methods, bilinear interpolation (BIL) method, bi-cubic interpolation (BIC) 
method, and curvature interpolation (CIM) method [9] have been compared. 
The least difference of PSNR between LSP method and CIM method is 0.38, 
the most is 3.18, and its average is 1.71. Comparing LSE with CIM and 
PSNR scale, the minimum difference is 0.13, the maximum is 3.12, and the 
average is 1.52. With regard to each row and comparing the results, it can be 
concluded that the proposed method of LSP has a better performance than 
other methods on the selected image. 
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Figure 6: Different images and standards for image magnification. 


Table 2 illustrates the results of comparing the suggested method and the 
other methods for the images in Figure 6, by SSIM criterion. The results show 
that the proposed method has a better performance than other methods with 
a good approximation. The minimum difference in SSIM for LSP by CIM is 
0.0011, the maximum is 0.02, and the average is 0.0096. In a comparison of 
LSE method with CIM and SSIM criterion, the minimum difference is 0.001, 
the maximum is 0.099, and the average is 0.0095. 

From a visual point of view, for example, the magnified image of Elaine 
is shown by LSP and LSE proposed methods, by BIL, BIC methods, as well 
as CIM method in Figure 7. Obviously, the images of the proposed method 
are more transparent and have less blur and good performance on the edges. 


Main image 


CIM LSP LSE 


Figure 7: The results of Elaine’s image magnification by the factor of 2 and using different 
methods. 
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Table 1: Comparison of the least square surface method and other methods by PSNR 


criterion. 


Methods | BIL BIC | CIM | Proposed method | Proposed method 
LSP LSE 

lena 30.20 | 30.11 | 30.58 33.47 33.34 
Elaine | 31.26 | 31.08 | 31.49 34.67 34.61 
Disk 30.94 | 30.82 | 31.56 33.07 32.70 
Fire 31.66 | 31.71 | 31.82 32.20 31.95 
Swan 31.19 | 30.74 | 31.44 33.22 32.73 
Dog 31.83 | 31.37 | 32.05 33.94 33.63 
Balloons | 31.29 | 31.73 | 32.28 33.52 34.07 
House 32.64 | 32.60 | 32.97 34.00 33.83 
Airplane | 29.70 | 29.59 | 30.12 31.34 31.30 
Fruits 26.17 | 26.08 | 26.49 28.44 28.20 
Average | 30.69 | 30.58 | 31.08 32.79 32.60 


Table 2: Comparison of the least square surface method and other methods by SSIM 


criterion. 


Methods BIL BIC CIM | Proposed method | Proposed method 
LSP LSE 

lena 0.9504 | 0.9469 | 0.9514 0.9714 0.9713 
Elaine | 0.9536 | 0.9479 | 0.9578 0.9589 0.9588 
Disk 0.9908 | 0.9908 | 0.9922 0.9956 0.9953 
Fire 0.9770 | 0.9767 | 0.9779 0.9857 0.9854 
Swan 0.9737 | 0.9727 | 0.9758 0.98420 0.9841 
Dog 0.9701 | 0.9691 | 0.9708 0.9814 0.9813 
Balloons | 0.9815 | 0.9808 | 0.9856 0.9893 0.9894 
House 0.9771 | 0.9763 | 0.9777 0.9867 0.9866 
Airplane | 0.9698 | 0.9696 | 0.9703 0.9831 0.9831 
Fruits 0.9518 | 0.9513 | 0.9523 0.9714 0.9712 
Average | 0.9696 | 0.9682 | 0.9712 0.9808 0.9807 
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5 Conclusion 


In this paper, a new method of the least squares surface has been used to 
enlarge images. Despite simple implementation and low computational com- 
plexity, this method provides more satisfactory results than bilinear interpo- 
lation methods, such as bicubic and curvature do. The study, has been con- 
ducted using one of the edge- improving techniques. The proposed method 
can be elaborated on in future studies through other edge- improving methods 
or nonlinear methods such as bivariate interpolation or radial basis functions 
interpolation with edge- shaving techniques. 
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